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ABSTRACT 

The recent discovery of a three-planet extrasolar system of HR 8799 by Marois et al. is a 
breakthrough in the field of the direct imaging. This great achievement raises questions on 
the formation and dynamical stability of the HR 8799 system, because Keplerian fits to as- 
trometric data are strongly unstable during ~ 0.2 Myr. We search for stable, self-consistent 
Af-body orbits with the so called GAMP method that incorporates stability constraints into the 
optimization algorithm. Our searches reveal only small regions of stable motions in the phase 
space of three-planet, coplanar configurations. Most likely, if the planetary masses are in 10- 
Jupiter-mass range, they may be stable only if the planets are involved in two- or three-body 
mean motion resonances (MMRs). We found that 80% systems found by GAMP that survived 
30 Myr backwards integrations, eventually become unstable after 100 Myr. It could mean that 
the HR 8799 system undergo a phase of planet-planet scattering. We test a hypothesis that the 
less certain detection of the innermost object is due to a blending effect. In such a case, two- 
planet best-fit systems are mostly stable, on quasi-circular orbits and close to the 5:2 MMR, 
resembling the Jupiter- Saturn pair. 
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1 INTRODUCTION 



The HR 8799 planetary system was detected by Marois et al. (2008) 
through the direct imaging. Soon, a new observation was added 
by Lafreniere et al. (2009) who reanalyzed images done in 1998, 
extending the observational window to ~ 10 years and four dif- 
ferent epochs. [We skip the most recent observation in (Fukagawa 
et al. 2009) that appeared after we finished this paper, because it did 
not change the initial condition]. Although the semi-major axes are 
large (about of 24, 36 and 68 au, respectively), the massive compan- 
ions strongly interact mutually. As we show, their orbits remain in 
extremely chaotic zone spanned by low-order MMRs. We attempt 
to constrain the initial conditions by available astrometric data and 
seemingly obvious requirement of astronomical stability. Our work 
complements papers of [Fabry cky & Murray^C lay (2008} and |Rei-| 
jdemeister et al.|p009) . Here, we follow a different approach that 
relies on quasi-global, self-consistent search for stable best-fit sys- 
tems, the so called GAMP (e.g., Gozdziewsk i et al.|2008| ) which 
was used to model the radial velocity data. The direct imaging 
seems also a particularly good target for this numerical technique. 

Following astrometric estimates of the semi-major axes, we 
see that the observational window covers a tiny part of orbital pe- 
riods which are counted in hundreds of years. The initial condition 
may be determined with a significant error. To illustrate this uncer- 
tainty, we map the multi-cube of astrometric coordinates x(t),y(t) 
and velocities v x (t),v y (t) [from the slope of x(t),y(t)] within 3a- 
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level of the linear model of x(t),y(t) onto osculating Keplerian ele- 
ments at the epoch of Sept. 18, 2008. This most reasonable choice 
follows the very short time-span of observations. The data set con- 
sists of 13 mean positions in [E, A/] -axes in ([Marois et al.|2008 > as 
well as one observation in ( Lafreniere et al. 2009|>; we also adopted 
a standard HIPPARCOS distance to the star of 39.4 ± 1 . 1 pc. 

The astrometric model is parameterized by the stellar mass 
mo, N tuples of Keplerian elements = (m[m Jup ], a[au], e, co[deg], 
^o[deg]), i.e., the mass, semi-major axis, eccentricity, argument of 
pericenter, and the mean anomaly (or longitude X), for each planet 
p= 1, . . . ,iV, respectively, and two Euler angles describing the in- 
clination (i) and the nodal longitude (£1) of the orbital plane with 
respect to the plane of the sky. The (Xv) 1 ^ 2 -function is build up 
from deviations of astrometric measurements from coplanar, pro- 
jected orbits. It depends indirectly on the astrophysical mass con- 
straints through the transformation of the velocity-Keplerian ele- 
ments. Following Marois et al. (2008), the planetary masses are: 
m b = 7^2 nij U p, m c = 10 ± 3 m Jup , m d = 10 ± 3 m Jup ; the mass of 
the parent star is mo = 1.5 ±0.3 hiq. They are roughly consis- 
tent with the recent, independent estimates of Reidemeister et al. 
( 2009). In our simulations, all masses are free parameters which 
are varied within their la-error ranges. (Moreover, a proper mass 
determination may be critically important for the stability analy- 
sis). Figure [l] shows levels of (Xv) m selected two-dimensional 
planes of osculating elements at the epoch of Sept. 18, 2008. The 
best-fit solution is marked with a green triangle. It is roughly con- 
sistent with a model of not too eccentric, face-on orbits by the dis- 
covery team. We found that the limited astrometric data permit a 
continuum of models with different orbital characteristics, e.g., ec- 
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Figure 1. Global topology of (Xv) 1 ^ 2 projected onto selected planes of Keplerian, osculating elements at the epoch of Sept. 18, 2008. Symbols are for 
the best- fits systems analyzed in this work. They are labeled accordingly with Table 1. Shaded areas are for la, 2a, and 3a-levels labeled in panels (a,d,f) 
[(Xv) 1 ^ 2 < 1-67, (X?) 1//2 < 1-85, and (x 2 ) 1 ^ 2 < 2.09, respectively] of the best Fit I marked with green triangle; darker shade means better fit. Black regions are 
for (x 2 ) 1 / 2 only marginally worse from the best- fit value. Other fits analyzed in this work are marked with white diamond (Fit II, unstable Trojans), yellow 
rectangle (Fit III, the best-fit stable configuration), and red circle (Fit IV, stable Trojans). See Fig.[2]for a geometry of these solutions. 



centricities within la-level of the best fit may be as large as ~ 0.4. 
The orbital periods consistent with relatively small (Xv) 1 ^ 2 ma Y ^ e 
found in a proximity of numerous low-order MMRs. In turn, these 
factors strongly affect the dynamical stability of the system. 

The geometry of the nominal, best-fit solution with (Xv) 1 ^ 2 ~ 
1.55, over-plotted on the original image, is illustrated in Fig.|2jl). 
This best fit system appears strongly unstable and self-disrupts af- 
ter ~ 0.2 Myr, so our conclusion is the same as in ( Fabry cky & 
Murray-Clay 2008). Moreover, according with our Fig. [T}l, very 
different orbital solutions are possible. For instance, two inner plan- 
ets might be involved in lc:ld MMR, or other low-order MMRs. 
An example of unusual Trojan configuration with only marginally 
worse (x 2 ) 1 ^ 2 ~ 1 -56 from the best-fit model is shown in Fig.^II). 
Also such "raw", kinematic fits are catastrophically unstable during 
the first Myr. Curiously, in these solutions, mo tends to the lowest 
possible limit that might indicate an internal inconsistency of the 
model with the data, if we recall that the stellar mass is constrained 
a priori. 

2 THE BEST-FIT STABLE CONFIGURATIONS 

Recent works (e.g.,|Juric & Tremaine|2008[|Chatterjee et al.|2008[ 
Scharf & Men ou|2009| and references therein) showed that com- 
pact planetary systems may evolve towards configurations span- 
ning wide ranges of orbital elements. The long-period planets could 
place constraints on early stage planet formation scenarios. An in- 
terpretation of the direct imaging surveys is also closely related to 
models of the dynamical relaxation ( Veras et al. 2009). Hence, even 
apparently odd solutions (like the Trojan configurations), consistent 
with observations, should not be skipped a priori. Because the par- 
ent star may be very young (30 Myr or less), and the dynamical 
separation of planets in terms of the mutual Hill radii, (Chatterjee 
|et al.|2008] >, K ~ 2, the three-planet system is strongly unstable in a 
few hundred orbital periods time-scale (see their Fig. 29; although 
these calculations are for more compact, Solar-system like models 
and planets in Jupiter mass range). So the HR 8799 system might be 
not yet dynamically relaxed, remaining in a stage of planet-planet 
scattering. On the other hand, we may be "fooled" by the significant 



errors of the initial condition implied by short time-base of the as- 
trometry. Then the requirement of the dynamical stability may help 
us to find long-living systems close to apparently unstable best-fit 
configurations. 

As is well known, the phase- space of a compact multi-planet 
system has non-continuous structure with respect to any notion 
of stability. The permitted region in the 18 -dimensional parame- 
ter space of the HR 8799 system is large and has complex shape. 
To explore it efficiently, we apply a variation of the GAMP method 
(see e.g., Gozdziewski et al. 2008 for details) which relies on self- 
adapting optimization based on the genetic algorithms (GAs) (e.g., 
Charbonneau 1995, Deb et al. 2002 ) and on "penalizing" unsta- 
ble configurations by an appropriate term added to the mathemat- 
ical value of (x 2 ) 1 ^ 2 - Here, the penalty term is expressed through 
the diffusion of fundamental fr equencies (Robutel & Laskar 2001 ; 
Sidlichovsky & Nesvorny 1996| >. 



An extensive GAMP search revealed long-term stable best- 
fit III (Table 1) illustrated in Fig.^III). We note that its (x?) 1/2 ~ 
1.88, remaining within 2a-range of the nominal, kinematic Fit I. 
To understand this solution, we computed its dynamical maps in 
terms of the Spectral Number (SN), the fast indicator invented by 
Michtchenko & Ferraz-Mello (2001), and the maxe (the maximal 
eccentricity attained during prescribed integration time). The SN 
map is illustrated in Fig[3] Fit III lies inside a small island of regu- 
lar motions (its width for the innermost planet is only ~ 0.3 au). A 
map of the max e indicator (not shown here) reveals that outside this 
region, one of orbits become highly eccentric that leads to catas- 
trophic events during ~ 1 Myr. Fit III describes a configuration in- 
volved in the Laplace-type three-body resonance, ld:2c:4b MMR. 
Its critical argument is shown in the middle panel of Fig. [4] A simi- 
lar solution was already found by Fabry cky & Murray-Clay (2008) 
and analyzed in more detail by Reidemeister et al.| ( |2009| >. Actu- 
ally, our Fit III is also unstable but on a very long time-scale. Af- 
ter ~ 400 Myr, the innermost eccentricity suddenly grows and the 
Laplace resonance disrupts (see two upper panels in Fig. [4), in- 
dicating a collision. Hence, the small amplitude of the resonance 
angle does not protect the system from the collision. In fact, Fit III 
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Figure 2. Orbital geometry of the best-fit configurations projected onto the plane of the sky and the true image of the system by combined photographs taken 
with the Keck II adaptive optics. The planets appear as red dots around the residual scattered light of the star (seen in the center). The best-fit osculating orbits 
at the epoch of Sept. 18, 2008 are drawn with yellow ellipses. Stright lines are for the apsidal lines of these orbits. The HR 8799 image credit: NRC Canada/C. 
Marois. See Table 1 for the osculating elements of these fits labeled accordingly. 
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Table 1. Osculating elements of the best-fit solutions at the epoch of Sept. 18, 2008. Formal errors of kinematic Fit I and II can be estimated graphically, 
see Fig.[T] Note that formal errors of dynamical GAMP Fits III and IV can be taken the same as for Fits I and II but the error ranges are strictly limited to 
stable regions of the phase space. The same concerns Fit V (see the text). The orbital inclination and nodal longitude (Fig.[T£) and arguments of pericenter 
(not shown here) are unconstrained. We set i G [0, 30°], consistent with estimates of (Marois et al.|2 008 ), following the rotation model of the star. 
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Figure 3. The SN map around stable Fit III (yellow rectangle in Fig.[TJ in 
the (a c ,e c ) -plane; other elements are fixed at their nominal values (see Ta- 
ble 1). Yellow colour is for strongly chaotic systems, black is for stable so- 
lutions. The maximal integration time for each pixel is 10 Myr (~ 25000Pd)- 



is formally chaotic that is indicated by the MEGNO in Fig. [4] This 
shows that the stability depends on long-term effects of the three- 
body mutual interactions and is tightly related to formally chaotic 
or regular character of tested configurations. 

Besides Fit III, we also found a stable fit related to 
lc:l d M MR, with moderate £ cd ~ 0.25 and a cd ~ 31 au, see 
IV). Its (Xv) 1 ^ 2 ~ 2 is still acceptably small because it lies 



Fig 

within the formal 3a-level of the best Fit I. The dynamical SN map 



of this fit is shown in Fig. [5] It reveals also a small island of stable 
motions having the width comparable to Fit III. Simultaneously, 
planet b remains in a narrow island close to (lc: Id): 3b MMR. The 
Trojans live at least over 3 Gyr — this system is close to quasi- 
periodic one, as indicated by (Y) (t) ~ 2 ( [Cincotta et al. 2003) over 
large part of the integration time (Fig. [6}. Still, a 10 Myr MEGNO 
map (not shown here) shows that the island of regular solutions 
is very tiny (~ 0.01 au). This fit is also weakly-chaotic although 
during first 600 Myr it appears as regular. This solution has pecu- 
liar small-amplitude librations of apsidal angle AU5(t) = G5 d — G5 C 
around 100° (the upper panel in Fig. [6). It might be the first case 
of asymmetric librations in the 1:1 MMR observed in a real sys- 
tem, and predicted already in low-order resonances (in particular, 
in 2:1 MMR, [Hadjidemetriou|2006] > . 

3 THE ASTRONOMICAL STABILITY OF THE SYSTEM 

The phase space of the HR 8799 system appears strongly chaotic, 
with tiny islands of regular two- and three-body MMRs. Hence, to 
study its long-term (but finite) evolution, we might rely on a notion 
of the astronomical stability ( |Lecar et al.|2001| ), rather than on the 
formal Arnold's stability analyzed above. The astronomical stabil- 
ity may be investigated only by the direct numerical integrations. 
Because the fate of chaotic configurations is hardly predictable, we 
attempted to gather statistics on initial conditions providing long- 
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Figure 7. Osculating elements of configurations astronomically stable over At = +100 Myr after the initial epoch to in terms of orbital periods ratio (panel 
a). Panels (b) and (c) illustrate the distribution in the (a, e) -plane of the final, dynamically relaxed systems of two planets. 
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Figure 4. The innermost eccentricity in Fit III (the top panel), the argu- 
ment of the Laplace resonance = — 3^ c + 2A,b (the middle panel) and 
MEGNO, (Y) (t) (the bottom panel). (Y) (t) converges to 2 for regular sys- 
tems and diverges linearly for chaotic motions ( Cincotta et al. 2003). 
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Figure 5. Dynamical SN map of a stable lc: Id MMR (Fit IV, see Table 1). 



living configurations, i.e., characterized by the event time Te of a 
close encounter/ejection of a planet from an initial system. 

We tested initial conditions within formal 3a-level of the nom- 
inal Fit I. To search for long-living systems, we again applied the 
GAMP algorithm, with the penalty term A(Xy) 1 / 2 multiplied by 
T = 1 — || Je/A^||, where At is the maximal integration time relative 
to the initial (present) epoch t = tQ. In the first simulation, we inte- 
grated the system over At = — 30 Myr (backwards), keeping track 
of solutions that survived as three-planet configurations. Next, the 
GAMP sample of ~ 2000 solutions "living in the past" was in- 
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Figure 6. Long-term stability of the best-fit solution with two inner planets 
involved in 1:1 MMR (Fit IV, see Tab. 1). The upper panel is for apsidal 
angle AG3 = G3d - G3 C . The bottom panel is for the MEGNO indicator. 



tegrated up to At = +100 Myr. The results are shown in Fig. [7] 
Panel [7^ is for the orbital periods ratio of systems that survived 
as three-planet configurations. It indicates that such systems are 
close to the Id: 2c MMR; the outer planets may be also involved 
in lc:2b MMR or other low-order MMRs (like 2c :5b). That agrees 
well with the results of Fabrycky & Murray-Clay ( 2008). Moreover, 
most of the tested systems self-disrupted. Two remaining panels in 
Fig. [7] are for the final osculating elements in the two-planet sam- 
ple. Due to intensive planet-planet scattering (we recall that K~2), 
the distribution spans large ranges of semi-major axes and almost 
whole available range of eccentricity. The strongly chaotic char- 
acter of the HR 8799 system leads to rapid collisions/ejections in 
most of tested configurations during at most a few Myr. That con- 
firms globally that the dynamical maps shown in Figs. |3|5| represent 
a generic picture of the phase space, although they were computed 
for particular (resonant) initial conditions. Still, only ~ 400 systems 
of the total population, i.e., less than 20%, survived the integrations. 
In fact, the sample is "biased" by the selection of systems surviving 
the integrations backwards. We found that the direct Monte-Carlo 
integrations leave much less than 1% of astronomically stable con- 
figurations after 100 Myr. Hence, the self-adapting GAMP is cru- 
cial in this test because the direct Monte-Carlo simulations would 
lead to unacceptable CPU overhead. Actually, there is no guarantee 
that systems astronomically stable in the 100 Myr test will also re- 
main stable on longer time-scale, as shows the case of Fit III. In this 
experiment, we also found ~ 20% of single-planet systems, and the 
rest in the sample ended as two-planet configurations. These dy- 
namically relaxed two-planet systems appear highly hierarchical, 
with a strong maximum of semi-major axes ratio a ~ 0.05 (see 
also Fig.[7]},c). 
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Figure 8. Dynamical (a c ,e c )-map of the best Fit V (see Table 1) of two- 
planet configuration in terms of the SN indicator. Its position is marked by a 
crossed circle. Most prominent structures of low-order MMRs are labeled. 



4 TWO-PLANET HYPOTHESIS 

Up to now, we assumed that the HR 8799 hosts three planets. Due 
to short observations that revealed planet d (a few weeks only), its 
orbit is unconstrained. Marois et al. (2008 ) claim a 6a detection 
of the common proper motion, consistent with the Keplerian orbit. 
It is enforced by the absence of HR 8799d in the HST images in 
1998 (Lafreniere et al. 2009 ) that otherwise should be seen in the 
substracted light annulus (C. Marois, private comm.). Still, we look 
here for an alternative explanation of the strongly unstable HR 8799 
system due to projected brown dwarf or already ejected planet that 
would be really too distant to influence orbits of HR 8799b,c. 

We repeated the GAMP experiment for such two-planet 
model. We found easily rigorously stable solutions with (Xv) 1 ^ 2 
comparable with the nominal, kinematic best-fit system. Elements 
of the best-fit solution are given in Table 1 (Fit V). The two-planet 
fits within la-bound span a wide range of semi-major axes ~ 10 au 
and may be stable up to e^ c ~ 0.15. Their orbits are initially close 
to anti-aligned ones. Planetary masses in such systems remain in 
the 10 m Jup range that is well consistent with astrophysical con- 
straints given in ( [Marois et al.||2008] >; we note that stable three- 
body fits tend to much lower masses than declared (Fabrycky & 
Murray-Clay 2008 ). Also the dynamical map in Fig. [8] shows ex- 
tended zones of stability and a proximity of the best-fit solution to 
the 5:2 MMR, recalling the Jupiter- Saturn pair in the Solar-system. 



5 CONCLUSIONS 

The dynamical analysis of available astrometric data of HR 8799 
reveal that its massive companions are involved in heavy mu- 
tual interactions. Assuming la-range of the planetary and star 
mass astrophysical estimates, the search for stable (regular) sys- 
tems brings only narrow and very limited islands of ordered mo- 
tions. Most likely, the system can be long-term stable if is in- 
volved in low-order two- or three-body MMRs (particularly, in 
the Laplace ld:2c:4b MMR). Here, we confirm the results of |Fab-| 
|rycky & Murray-Clay | ( [2008] ) and |Reidemeister et al.| ( [2009| >, which 
we derived after independent, quasi-global GAMP calculations. 
Moreover, also peculiar ld:lc MMR Trojan systems stable over a 
few Gyr can be found. 

The outstanding discovery, in the light of the dynamical anal- 
ysis, brings a few open questions. How the three-planet system may 
be captured in such tiny regions of stable motions? Are in fact plan- 
etary masses much lower than estimated? Or is the system sub- 
stantially non-coplanar? Both these factors could extend the zones 
of stability. While the masses may be constrained by astrophysical 



factors and astrophysical-age estimates (Marois et al. 2008 ; Reide- 
meister et al. 2009 ), we can say little on the real mutual inclinations. 
Further, if we "skip" the less constrained object, the sub-system of 
outermost planets is stable, resembling the Jupiter- Saturn pair, even 
if the masses are large, apparently solving the puzzle. It may be ver- 
ified soon, thanks to the shortest orbital period of planet d. 

Actually, should we expect that the system is or must be sta- 
ble? Its parent star is very young, and we may have an opportu- 
nity to observe a system undergoing the dynamical relaxation. The 
statistical analysis suggest, that the final fate of coplanar systems 
constrained by available astrometric data most likely will be two- 
planet, highly hierarchical configuration with eccentric orbits. Our 
calculations show that less than 20% of systems stable in the past 
and remaining in the neighborhood of the best stable Fit III remain 
stable after 100 Myr. Likely, even much less number of configu- 
rations survive longer time due to possible, chaotic effects of the 
three-body interactions [MMRs overlapping, (Murray & H olman| 
|2001| )]. A conclusion of |Fabrycky & Murray-Clay (2008 ) may be 
repeated here. Although the HR 8799 has been directly imaged, 
the interpretation of its images is very difficult and yet non-unique. 
Longer observations are required to constrain orbits of its planets. 
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